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Abstract 



A method for moving least squares interpolation and differentiation is 
presented in the framework of orthogonal polynomials on discrete points. 
This yields a robust and efficient method which can avoid singularities 
and breakdowns in the moving least squares method caused by particular 
configurations of nodes in the system. The method is tested by applying 
it to the estimation of first and second derivatives of test functions on 
random point distributions in two and three dimensions and by examining 
in detail the evaluation of second derivatives on one selected configuration. 
The accuracy and convergence of the method are examined with respect 
to length scale (point separation) and the number of points used. The 
method is found to be robust, accurate and convergent. 

1 Introduction 

The moving least squares method [3, chapter 7] is a technique for interpola- 
tion [6] and differentiation [1,2,7-10,13] on scattered data. The purpose of 
this paper is to examine the moving least squares problem in the framework of 
orthogonal polynomials, as applied to the estimation of derivatives. 

In applications of the type considered here, the data supplied are the posi- 
tions of N points Xi, i — I, . . . , N and corresponding values /j. At one of these 
points the derivative is to be estimated. This is done using an interpolating 
polynomial P(x) which minimizes the error: 



where Wi is a strictly positive weight. The polynomial P(x) can be computed 
by direct solution of a least squares problem and then used to interpolate /(x) 
or to estimate its derivatives. 
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There arc a number of applications where moving least squares is used to 
estimate derivatives of a function specified at discrete points. One is the esti- 
mation of gradients of vorticity in Lagrangian vortex methods [7, 8] , where the 
gradients are estimated in two or three dimensions by fitting a second order 
polynomial to the components of vorticity and differentiating the polynomial. 
It was noted that "when computational points become very isolated, due to 
inadequate spatial resolution, the condition number of the matrix [used in fit- 
ting the polynomial] becomes very large" [7]. The solution proposed for this 
ill-conditioning was to add additional points to the fit. It appears that this 
problem may have been caused by another effect which has been noted by au- 
thors who use moving least squares to solve partial differential equations on 
irregular meshes or using mesh-free techniques. 

In the work of Schonaucr and Adolph [9,10], a finite difference stencil is 
developed using polynomials which interpolate data on points of an unstructured 
mesh. The points used in the polynomials are selected by choosing more points 
than there are coefficients in the fitting polynomial because "in m nodes usually 
there is not sufficient information for the m coefficients" [9] or, restated, "there 
are linear dependencies on straight lines" [10] . The number of extra points used 
in fitting the polynomial was determined through experience and testing. This 
raises the issue of the arrangement of the points used in deriving a polynomial 
fit. 

The issue has been addressed recently by Chcnoweth et al. [2] who considered 
the problem of how to find a least squares fit on points of an unstructured mesh 
in order to generate a stencil, while avoiding singularities caused by particular 
point configurations, a general form of the problem of "linear dependencies" [10]. 
They state the conditions under which such singularities will arise and state a 
criterion determining when it will not be possible to make a least squares fit of 
a given order on a given set of points in two dimensions. This will happen when 
selected points are spanned by the same polynomial, for example, when fitting 
a second order polynomial to points which lie on an ellipse in two dimensions. 
They also give an algorithm for a moving least squares fit which determines 
when more points must be added in order to avoid singularities, and which 
additional points will be useful. 

Another recent paper employing moving least squares methods for three- 
dimensional meshless methods [13] proposes an approach which may help avoid 
the problem of singularities. The method is to derive a set of basis functions 
which are orthogonal with respect to an inner product defined on the set of 
points. The use of orthogonal functions has the advantage of improving the 
condition number of the system to be solved to form the least squares fit and, 
in this case, allows a smaller number of basis functions to be used. The authors 
do not, however, discuss the problem of singular point configurations other than 
stating the number of points included in the fit must be large enough to make 
the system matrix regular, which corresponds to the avoidance of singular or 
ill-conditioned arrangements of nodes. 

Strangely, there does not yet appear to be a published moving least squares 
method which explicitly frames the problem in terms of orthogonal polynomi- 
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als. The aim of this paper is to present a method using resuhs from the theory 
of orthogonal polynomials in multiple variables [11] to restate the problem in 
a manner which detects singular point configurations and generates a set of 
orthogonal polynomials which arc unique for the points considered. The poly- 
nomials derived can then be used directly in computing a fit to the function on 
the specified data points. The method is quite general and does not require a 
knowledge of which configurations give rise to singularities. In three or more 
dimensions these configurations are not easily visualized and, furthermore, a 
singular value decomposition becomes increasingly expensive. 

2 Discrete orthogonal polynomials for scattered 
data 

The theory of classical orthogonal polynomials of several variables is well-developed [11] 

but that of polynomials orthogonal on discrete points is not as advanced. A 
recent paper [12], however, establishes basic properties of discrete orthogonal 
polynomials and gives algorithms for their derivation. In particular, it estab- 
lishes the theoretical foundations which allow us to say, given a set of points, 
whether orthogonal polynomials of a given order exist on these points and, if 
they do, what those polynomials are. In this section, we will summarize the 
mathematical tools required to derive and apply polynomials orthogonal on dis- 
crete points. We use the standard notation in which a polynomial of several 
variables is defined as a weighted sum of monomials: 

n 

P(x)=^^,x-, (2) 

where x = {xi,X2, ■ ■ ■ jXa), x e M'', ot = {ai,a2, ■ ■ ■ ,ad), a G Ng and the 
monomial terms x" = 11^=1 ■ The degree of P(x) is maxjajj where jaj = 

2.1 Generation of orthogonal polynomials 

The first basic tool required is a scheme to generate a set of polynomials which 
are orthogonal on a given set of points with respect to some weight function. 
This can be done using standard matrix operations [4, 5] using the procedure 
given by Xu [12]. First, we define the inner product: 

N 
i=l 

where / and g are functions evaluated at the data points Xj and Wi is the weight 
corresponding to Xj, with Wi > 0. 

The first step in generating the orthogonal polynomials is to find a set of 
monomial powers cxj which spans the polynomial space on Xj. This is done by 
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starting with the monomial 1 and systcmaticaUy adding monomials of increasing 
degree cXj. As each monomial is added, a matrix 



X = 



^1 ^2 • • • -^N 

J^Q ... JS~\T 



is generated for some initial value n. New rows are added to X for succes- 
sive values of otj, taken in lexicographical order at each \a\. The rank of X 
is checked at each step; if it is rank-deficient, the newly added monomial is 
rejected. Otherwise, it is added to the list of otj to be included in generating 
the polynomials. Rejection of a monomial power will happen because the point 
configuration is singular for the combination of monomials which would result 
from including the new ocj. Monomials are added until X is square and of full 
rank. The output of this procedure is a list of monomial powers which together 
span the polynomial space on the data points. 

To generate the orthogonal polynomials from the list of monomials, the 
following procedure is used: 

1. generate the symmetric, positive definite matrix M, with Mij = (x°' , x"^ ) . 

2. perform the decomposition M = SDS^, where D is a diagonal matrix 
and S is lower triangular. 

3. solve S'^R = D-i/^ where D-^/^ = dia.g{{diwi)-^/^, . . . ,{dNWN)-^/^}. 
This can be done using an LU solver with rearrangement of the matrix 
entries. 

4. the matrix R now contains the coefficients of the orthogonal polynomials. 

In implementing the method, we note that S"^ can be found directly by using 
the algorithm given by Golub and Van Loan [5, page 138] with the row and 
column indices switched. 

The orthogonal polynomials Pi are now: 

n 

Pi(x) = ^i?,,x«^ 

and for later convenience, we scale the coefficients on the inner products (Pi(x)Pj(x)) 
to give an orthonormal basis. 

2.2 Fitting data on sets of scattered points 

Given the set of orthogonal polynomials -Pi(x), generation of a least-squares fit 
is trivial. By orthogonality: 



/(x)«.^c,P,(x), (3) 
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where the constants Cj are given by: 

3 

Rearranging to give the interpolant as a weighted sum over the data points: 

/(x)«5^^;,/„ (4) 

3 



=WjY,Pi{^3)Pi{^)- 



Derivatives of /(x) can also be estimated as a weighted sum of the function 
values at the points of the distribution, to generate a differentiation stencil: 



3 

Ql+m+--- 



(Ir, 



with the derivatives of Pi(x) being computed directly from the coefficients in 
the matrix i? of §2.1. 

In summary, a derivative of a function /(x) given on a set of points can be 
estimated at some point xq using these steps: 

1. select N points in the region of xq, including Xq itself; 

2. generate a set of orthonormal polynomials for the selected points, using 
the procedure of §2.1; 

3. evaluate the weights fj'™ given by equation 5; 

4. calculate the derivative as the weighted sum of the function values. 

An important point is that strictly this procedure can only evaluate linear 
combinations of derivatives. In an extreme case, where only two points are used, 
the available monomials will be 1 and Xi (or X2 depending on the lexicographical 
ordering used). This allows linear interpolation of a function / between the two 
points and estimation of the derivative on a straight line joining them. This 
derivative will be a linear combination of df /dxi and df /dx2, with the precise 
combination depending on the orientation of the two points. In practice, this 
should not be a serious limitation since the monomials used in the polynomials 
are known and it is possible to determine whether there is a full set available 
for the determination of all derivatives of a given order. 
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3 Performance 

To illustrate the operation of the method, the first results presented are orthog- 
onal polynomials on regular arrangements of points. The first is a regular 3x3 
grid in (— f , 1) x (—1, 1). Upon application of the algorithm of §2.1, the mono- 
mials which form the final matrix X are 1, xi, X2, xf, X1X2, x^, x\x2-, X1X2 and 
a;f a;| and the resulting orthogonal polynomials are: 

Po = l, (6a) 

^1 = 21/231/2 ^1' (6b) 

= ^IT^V^^^, (6c) 
21/2 I 

^3 = -— + (6d) 

P4 = -X1X2, (6e) 



2 

2^/2 1 „ 
^5 = -— + ^4 (6f) 

1 31/2 

^6 = "^'^2 + —xlx2, (6g) 
1 3I/2 

P'^ = + —xixl, (6h) 

P8 = ^-xl-xl + ^xlx2- (6i) 

If the procedure is applied to six points equally spaced on a unit circle, the 
resulting polynomials are: 

^0 = 2uk^' (^^) 

^1 = ^^1' (7b) 
P2 = ^X2, (7c) 

^3 = -^ + ^^?> (7d) 
2 

Pi = -^XiX2, (7e) 
3I/2 21/2 
= -^^1 + 2^^?- (7f) 

A number of general issues are illustrated by these examples. The first is 
the obvious one that there are no more polynomials than there are points. This 
means that although the polynomials are notionally up to third order in both 
cases, in practice neither group of functions has a complete set of monomials 
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capable of spanning all polynomials up to cubic. Secondly, if the orthogonal 
polynomials can be generated, there is no benefit in adding extra points once 
a complete set of functions is available: the result of adding more points is to 
yield an incomplete set of higher order polynomials. In applications, it may well 
be better to have a lower order, but complete, system to fit the functions on the 
points. 

3.1 Random point distributions on the unit disc or ball 

The first set of results presented are average data for a large number of tests 
conducted with varying order and length scale. Following the example of Bel- 
ward et al. [1], the accuracy and robustness of the computational scheme are 
tested by estimating the derivatives of a prescribed function using a set of points 
randomly distributed over the unit disc or ball. The functions used are: 



The functions have been chosen to give a function which can be fitted exactly 
by a polynomial (/i), a Gaussian of the type found in various applications such 
as vortex dynamics (/2) and a Gaussian weighted to give an asymmetry with 
a consequent non-zero derivative at the evaluation point (/a). The evaluation 
point was fixed at and N = 8, 16, 32, 64, 128 random points were distributed 
uniformly in angle and radius over the unit disc or ball. Unit weights = 1 were 
used. Given the values of /(x^), the algorithm of §2.2 was used to estimate first 
and second derivatives of the function at the evaluation point. To examine the 
convergence rate, the procedure was repeated by using the values f{axi), where 
cr, < (T < 1, is a scaling factor which has the effect of contracting the point 
distribution around the evaluation point. It is assumed that in applications, a 
point distribution will be scaled to a normalized radius and the result of the 
function evaluation rescaled afterwards, a procedure which is modelled using the 
scaling factor a. Tests were repeated on 32 different random point distributions 
and mean absolute errors estimated. The presented results are mean absolute 
error, mean number of rejected monomials and convergence rate with cr, for 
different values of N and different functions /(x in two and three dimensions. 

Two sample sets of results are shown graphically to illustrate the perfor- 
mance of the method, with the relevant performance parameters for all tests 
summarized in tabular form. 

Figure 1 shows the performance data for evaluation of dfi/dxi in two di- 
mensions. This function can be fitted exactly by a polynomial of sufficiently 
high order as is clear from the results. The first plot of mean absolute error 
against a, shows that all three orders give identical results. This is because, with 
only eight points available in the fit, the three systems of functions are identical. 



/i(x) 

/2(X) 
,/3(x) 




(8a) 
(8b) 
(8c) 



XiC 
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log2 a N 

Figure 1: Evaluation of df\/dx\ in two dimensions: mean absolute error against 
scale factor u for, reading left to right, N = 8, 16, 32, 64, 128 points; final plot 
mean number of rejected monomials Nr against number of points N. Second 
order fit shown solid, third order dashed, fourth order long dashed. 
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Table 1: Convergence rates for first derivatives in two-dimensional problems 







8 


16 


32 


64 


128 


(128) 
mill 




(128) 




/l 


2 


4.00 


4.00 


4.00 


4.00 


4.00 


5.21x10" 


10 


3.41x10 


-5 




3 


4.00 


4.00 


4.00 


4.00 


4.00 


3.66x10 


-8 


2.40x10" 


-3 




4 


4.00 


4.00 


4.00 


4.00 


4.00 


3.98x10- 


17 


2.61x10" 


12 




2 


2.00 


2.03 


4.00 


4.00 


4.00 


1.09x10 


-9 


7.18x10- 


-5 




3 


2.00 


2.05 


4.00 


4.00 


4.00 


3.03x10 


-9 


1.98x10" 


-4 




4 


2.00 


2.05 


4.00 


5.42 


5.42 


4.99x10" 


13 


1.24x10 


-6 


h 


2 


3.00 


3.00 


3.00 


3.00 


3.00 


2.18x10 


-7 


8.93x10 


-4 




3 


3.00 


3.00 


3.02 


4.80 


4.80 


1.96x10" 


11 


1.05x10" 


-5 




4 


3.00 


3.00 


3.02 


4.35 


4.35 


4.99x10" 


11 


6.22x10 


-6 



Increasing the number of points to sixteen, the second order fit is slightly better 
than the other two which are themselves identical. The final plot, showing the 
mean number of monomials rejected for each N explains why. With N = 8, all 
three fits are underspecified while at TV = 16, the second order fit has a full set 
of six monomials but the third and fourth order fits are still short of useable 
terms. As N increases to 64, the third and fourth order fits gain a complete set 
of monomials (ten and fifteen, respectively) and the full accuracy of the exact 
fourth order fit becomes clear. The third order fit's error, however, is larger than 
that of the second order, probably because the third order monomials, which 
do not fit the symmetric function /i, introduce spurious terms in the fit. The 
results for A'' = 64 and N = 128 are identical, because at this stage a full set 
of monomials is available for all of the orders considered and adding additional 
points gives no extra benefit. In all of the cases considered, the error reduces 
with a at the same rate, though with quite different error magnitudes, as will 
be seen in the tabular data presented later. 

Figure 2 shows the performance of evaluation of df2/dxf in three dimensions. 
For reference, there are 10, 20 and 35 monomials in a fully-specified polynomial 
of order two, three and four respectively. As in the two-dimensional case, the 
first two plots show identical behavior of the error for all three fits, due to the 
number of points being insufficient to generate a fully specified polynomial. The 
third plot, N = 32, shows the second order fit being slightly better than the 
other two, which are identical, since this is now fully specified. As the number 
of points is increased the three fits begin to separate, although none has a clear 
advantage over the others until N = 128, where the fourth order fit has a full 
set of monomials available, as shown in the final plot. The convergence rate, 
however, reduces at small a which may be due to floating point errors. In 
the fourth order fit the monomials are 0(x*) and the resulting inner products 
0{x^). When a = 2"^, the maximum value of .x* is 2"^^ w 2 x 10""'^'', for points 
furthest from the evaluation position: the corresponding term for those points 
nearest the evaluation position will be much smaller, comparable to the floating 
point precision of the computer. 

Numerical data summarizing the results of all of the tests carried out are 
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log2 o- N 

Figure 2: Evaluation of f2/dx\ in three dimensions: mean absolute error |e| 
against scale factor a for, reading left to right, = 8, 16, 32, 64, 128 points; 
final plot mean number of rejected monomials A^,. against number of points N . 
Second order fit shown solid, third order dashed, fourth order long dashed. 
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Table 2: Convergence rates for second derivatives in two-dimensional problems 







8 


16 


32 


64 


128 


(128) 
mill 




(128) 
f'-max 




A 


2 


4.00 


4.00 


4.00 


4.00 


4.00 


7.39 X 10 


-8 


4.84 X 10" 


-3 




3 


4.00 


4.00 


4.00 


4.00 


4.00 


1.34 X 10 


-7 


8.84 X 10- 


-3 




4 


4.00 


4.00 


4.00 


4.00 


4.00 


2.65 X IQ- 


15 


1.74 X 10- 


10 




2 


2.00 


2.02 


4.00 


4.00 


4.00 


6.25 X 10 


-8 


4.09 X 10 


-3 




3 


2.00 


2.03 


4.00 


4.00 


4.00 


4.10 X 10 


-6 


2.67 X 10- 


-1 




4 


2.00 


2.03 


4.00 


5.01 


5.01 


7.49 X 10- 


11 


6.02 X 10- 


-5 


h 


2 


3.00 


3.00 


3.00 


3.00 


3.00 


3.15 X 10 


-5 


1.29 X 10 


-1 




3 


3.00 


3.00 


3.00 


4.88 


4.88 


1.37 X 10" 


10 


9.64 X 10 


-5 




4 


3.00 


3.00 


3.01 


4.75 


4.75 


3.29 X 10 


-9 


1.52 X 10 


-3 



presented in Tables 1 4. In each tabic, the convergence rate of the fits is pre- 
sented for each of the three test functions at each value of N, the size of the 
point set. In addition, to compare the errors proper, the final two columns give 
the maximum and minimum mean errors for fits performed using 128 points. 
Convergence rates r were found by a least squares fit |e| = coa^. 

Table 1 shows the performance data for evaluation of dfi/dx\. The results 
are much as might be expected, with a small minimum error in each case, 
especially for /i, the polynomial and with smooth convergence for most fits 
on most functions. The exceptions are the fourth order fit to A and /s and 
the third order fit to f^. In these cases, at large point number > 64, the 
errors are small, as would be expected, but the convergence is not as smooth as 
expected. For /s the mean convergence rate of the fourth order fit is also less 
than that of the third order, although the absolute errors are comparable. This 
is probably due to a combination of the floating point issue mentioned earlier 
and the inability of the fourth order monomials to capture the behaviour of the 
function near 0. 

Table 2 shows the equivalent results for evaluation of d"^ fidx\. The results 
show the same trends as in Table 1, with the fourth order fit giving very small 
minimum errors for all three functions but with the third order fit being slightly 
superior for f^. 

Tables 3 and 4 give data for first and second derivative evaluation in three 

dimensions. The general trends are similar to those in two dimensions but there 
some differences worth noting. In evaluating the derivatives of /s, the fourth 
order fit, as in the two-dimensional case, does not perform as well as the third 
order, although the error is still small. The difference is in the convergence rate 
as N increases from 32. At = 64, the convergence rate drops from 3 to 2 
before increasing again to 4.39, in contrast to the behavior of the third order 
fit. This is probably because at iV = 32, the third and fourth order fits are both 
underspecified (see the final plot of Figure 2) but the third order fit gains a full 
set of monomials at A'' = 64. The fourth order fit has a full set of third order 
monomials but has rejected, on average, eight monomials, leaving only seven 
fourth order terms in the polynomials. This leads to error behavior which is 
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Table 3: Convergence rates for first derivatives in three-dimensional problems 







8 


16 


32 


64 


128 


(128) 
min 




(128) 
cmax 




fl 


2 


4.00 


4.00 


4.00 


4.00 


4.00 


3.14 X 10 


-7 


2.05 X 10" 


-'2 




3 


4.00 


4.00 


4.00 


4.00 


4.00 


1.73 X 10" 


-8 


1.13 X 10" 


-3 




4 


4.00 


4.00 


4.00 


4.00 


4.00 


4.55 X 10" 


14 


2.98 X 10" 


-9 


/2 


2 


2.00 


2.00 


4.00 


4.00 


4.00 


4.73 X 10 


-9 


3.09 X 10- 


-4 




3 


2.00 


2.00 


3.97 


4.00 


4.00 


1.18 X 10- 


-6 


7.79 X 10- 


-2 




4 


2.00 


2.00 


3.97 


3.99 


3.61 


9.13 X 10" 


10 


2.31 X 10" 


-5 


fs 


2 


3.00 


3.00 


3.00 


3.00 


3.00 


5.12 X 10 


-7 


2.09 X 10" 


-3 




3 


3.00 


3.00 


3.00 


4.75 


4.75 


1.50 X 10" 


10 


6.89 X 10- 


-5 




4 


3.00 


3.00 


3.00 


2.01 


4.39 


7.88 X 10- 


-9 


1.16 X 10" 


-3 



Table 4: Convergence rates for second derivatives in three-dimensional problems 







8 


16 


32 


64 


128 


(128) 
min 




(128) 
tmax 




fl 


2 


4.00 


4.00 


4.00 


4.00 


4.00 


1.66 X 10 


-a 


1.09 X 10" 


-1 




3 


4.00 


4.00 


4.00 


4.00 


4.00 


1.35 X 10" 


-6 


8.91 X 10" 


-2 




4 


4.00 


4.00 


4.00 


4.00 


4.00 


1.11 X 10- 


11 


7.32 X 10" 


-7 


/2 


2 


2.00 


2.01 


4.00 


4.00 


4.00 


2.32 X 10- 


-6 


1.51 X 10" 


-1 




3 


2.00 


2.01 


4.00 


3.98 


3.98 


6.25 X 10" 


-7 


3.85 X 10- 


-2 




4 


2.00 


2.01 


4.00 


3.85 


3.83 


9.17 X 10- 


-8 


3.95 X 10- 


-3 


/3 


2 


3.00 


3.00 


3.00 


3.00 


3.00 


7.14 X 10- 


-5 


2.91 X 10- 


-1 




3 


3.00 


3.00 


3.00 


5.00 


5.00 


4.23 X 10 


-9 


4.44 X 10- 


-3 




4 


3.00 


3.00 


3.00 


4.13 


3.50 


3.36 X 10 


-6 


4.56 X 10- 


-2 
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Table 5: Fourth order monomials included in polynomial fits on point set of 
Figure 3 and resulting error in d"^ j j dxidx^ 

N Fourth order monomials included log2 |e| 



worse than that from a fully specified third order polynomial. As N increases 

again, to 128, the fourth order polynomial is fully-defined and its convergence 
rate recovers, with the caveat that it may now be affected by floating point 
errors. 

To examine the behavior of the method in more detail, we look at the results 
of a test on a fixed random point configuration. The task is to estimate the 
second derivatives of /2 = exp(— i?^) in two dimensions with a = 1/2^ using 
a nominally fourth order fit. Twenty randomly distributed points, sorted by 
distance from the evaluation point x = 0, are used and a fit is generated using 
A'^ = 6, 7, . . . , 20 of these points in turn. The results are shown in Figure 3: the 
first plot shows the point configuration with points numbered by distance from 
the evaluation point, the second plot shows the number of monomials rejected 
at each N and the third shows the error in the estimate of each of the three 
second-order derivatives. 

The error behavior demonstrates some of the detailed features of using the 
estimation scheme. For a fully specified fourth order polynomial, fifteen mono- 
mials are required. It is only at = 20 that these all become available, with 
sufficient points being used to avoid singular configurations. For N = 6, 7, 
the error is quite large, due to the fit being a set of defective second order 
polynomials — ten monomials are rejected and only five terms are available; for 
a notionally second order fit. At VV = 8, one more monomial bc!comc;s available 
and the error drops immediately since this the method is now a fully specified 
second order fit to a set of points closc^ to the evaluation point. The error re- 
mains constant up to A'^ = 18, even though the number of points is increasing, 
since the additional terms are third order and do not contribute to a fit on the 
symmetric function in question. There is a drop in the error at = 18, fol- 
lowed by an increase and another reduction. This can be explained by looking 
at which monomials have been included or rejected in the fits. 

Table 5 shows the fourth order terms which are included in the polynomial 
fits whose error behavior is shown in Figure 3, for 17 < N < 20, with the final 
column showing the error in a second order derivative. Each of the fits has a 
full set of lower order monomials and, in principle increasing A'^ allows more 
fourth order terms to be added. In practice, as can be seen, on this point set, 
at N = 17, two fourth order terms have been rejected and the accuracy suffers. 
At A'' = 18, the term added, and the accuracy improves markedly: this 

monomial is symmetric and is useful in a fit to the symmetric function being 
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Figure 3: Single point configuration test: top: point configuration; middle: 
number of rejected monomials against number of points used; bottom: error in 
d'^f/dxl (solid), d'^f/dxidx2 (dashed) and d'^f/dxl (long dashed). 
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differentiated. At A'' = 19, the monomial x^x"^ is rejected and x}.T2 is included. 
This term increases the error as in seen in the final column of Table 5 and in 
Figure 3. A full set of monomials is only available for this point configuration 
whc;n A^ = 20, resulting in the full expected accuracy. 

This behavior is slightly unexpected given that the orthogonal polynomi- 
als derived for any value of N span the polynomial function space on those 
points. This raises the issue of which basis functions should be used in applica- 
tions. Chenoweth et al. [2] discuss the problem of singular point configurations 
in the context of the singular value decomposition of a matrix containing the 
monomials evaluated at each point. As these authors note, a singular value de- 
composition shows which basis functions span the null space of polynomials on 
the points, allowing the detection of singular point configurations. The opposite 
fact is also true: the singular value decomposition yields a set of basis functions 
which span the function space on the points and, indeed, will indicate which 
of these basis functions are best determined. The problem, as we see above, 
is that even when a full set of well-determined basis functions is available, it 
is not guaranteed that they form a suitable basis for the evaluation of deriva- 
tives, unless some extra measures are taken, as in Chenoweth et al's work where 
derivatives are included in the general form of the function to be fitted [2] . 

Given that in many applications, it will not be known in advance which 
terms will be most useful in fitting a function, it is recommended that only fully 
specified polynomials be employed, with the order depending on the accuracy 
required and the point density available. 

4 Conclusions 

A method for moving least squares interpolation and differentiation using or- 
thogonal polynomials has been presented and tested on random point distribu- 
tions. The method makes \isc of the theory of discrete orthogonal polynomi- 
als in multiple variables and deals with the problems caused by singular point 
configurations by adjusting the terms of the polynomial. It is concluded that 
the method is robust and capable of detecting and compensating for singular 
configurations. In applications, it is recommended that the highest order of 
polynomial for which a full set of monomials is available be used in computing 
derivatives. 

Acknowledgements 

The author thanks the anonymous referees who read the original submission 
with great care, making many useful comments on the method of the paper and 
on its presentation. 



Moving least squares via orthogonal polynomials 



16 



References 

[1] J. A. Belward, I. W. Turner, and M. Ilic, On derivative estimation 
and the solution of least squares problems, Journal of Computational and 
Applied Mathematics, 222 (2008), pp. 511-523. 

[2] S. K. M. Chenoweth, J. Soria, and A. Ooi, A singularity- avoiding 
moving least squares scheme for two-dimensional unstructured meshes, 
Journal of Computational Physics, 228 (2009), pp. 5592-5619. 

[3] G. E. Fasshauer, Multivariate meshfree approximation. Course notes, 
http : //www . math . iit . edu/~f ass/603 .html, 2003. 

[4] M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, 
M. Booth, and F. Rossi, GNU Scientific Library Reference Manual, 
Network Theory Ltd, Bristol, United Kingdom, 2005. 

[5] G. H. GoLUB and C. F. Van Loan, Matrix computations, Johns Hopkins, 
Baltimore, Maryland, third ed., 1996. 

[6] D. Levin, The approxim,ation power of moving least-squares. Mathematics 
of Computation, 67 (1998), pp. 1517-1531. 

[7] J. S. Marshall, J. R. Grant, A. A. Gossler, and S. A. Huyer, 
Vorticity transport on a Lagrangian tetrahedral mesh, Journal of Compu- 
tational Physics, 161 (2000), pp. 85-113. 

[8] C. Moussa and M. J. Carley, A Lagrangian vortex method for un- 
bounded flows, International Journal for Numerical Methods in Fluids, 52 
(2008), pp. 161-181. 

[9] W. SCHONAUER AND T. Adolph, How WE solvc PDEs, Journal of Com- 
putational and AppUed Mathematics, 131 (2001), pp. 473-492. 

[10] , Higher order may be better or may not be better: Investigations 

with the FDEM (Finite Difference Element Method), Journal of Scientific 
Computing, 17 (2002), pp. 221-229. 

[11] Y. Xu, Lecture notes on orthogonal polynomials of several variables, in 
Inzell Lectures on Orthogonal Polynomials, W. zu Castell, F. Filbir, and 
B. Forster, eds., vol. 2 of Advances in the theory of special fimctions and 
orthogonal polynomials. Nova Science Publishers, 2004, pp. 135-188. 

[12] , On discrete orthogonal polynomials of several variables. Advances in 

Applied Mathematics, 33 (2004), pp. 615 632. 

[13] Z. Zhang, P. Zhao, and K. M. Liew, Analyzing three-dimensional po- 
tential problems with the improved element-free Galerkin method. Compu- 
tational Mechanics, 44 (2009), pp. 273-284. 



